In this practical, we will visualize a few data sets with R.
Geobr library - we will open and visualize the Brazilian territory data using the library called Geobr. We will also visualize the Amazon region and all areas of indigenous land, conservation units and urban areas of Brazil.
Prodes data - we will open and visualize the deforestation of the Brazillian Amazon from 2000 to 2020. Deforestation of the Amazon or the other ecosystems of Brazil is freely available from the PRODES portal http://terrabrasilis.dpi.inpe.br/en/download-2/
Terrabrasilis - we will use Terrabrasilis Analytics API to access the data of the Terrabrasilis portal https://github.com/terrabrasilis/terrabrasilisAnalyticsAPI. TerraBrasilis contains an open-source map server (GeoServer (http://terrabrasilis.dpi.inpe.br/geoserver/web/) for displaying interoperable spatial data.TerraBrasilis also integrates customized mapping API https://github.com/Terrabrasilis/terrabrasilis-api procedures and utilities (https://github.com/Terrabrasilis/terrabrasilis-util) as a simple means of creating interactive maps. We will use Terrabrasilis to calculate the area of deforestation for states, for example, Para or Roraima states.
Local case study - we will evaluate the effects of forest fires and droughts from 2016 in a central area of the Amazon called Tapajos National Forest (TNF). First, we will open and visualize the deforestation around this area from different databases (Prodes and MapBiomas). Second, we will visualize the area burned in 2016 and the area of secondary forest still available in the region. Third, we will plot a few maps of precipitation data around this region to investigate the effect of droughts in this region. Fourth, we will overlay all roads of the region inside the land-use map of TNF. Fifth, we will divide the entire region with a few grids of 10 by 10 km (cells). With those cells, we are going to calculate the fraction of burned forest by cells, the number of active fires or hot spots of fires per cell, and the amount of accumulated water deficit. Sixth, with all these grid cells, we are going calculate the correlation between proportions of agriculture (ag_p), pasture (pa_p), urban area (ua_p), non-forest (nf_p), water bodies (wt_p), active fire from Jun 2015 to 2016 (af_total), and burn scar from Jun 2015 to Aug 2016 (bs_total). Finally, we will calculate the correlation between active fire from Jun 2015 to 2016 (af_total) and environmental factors that may have contributed to the droughts, such as drought season length (ds), minimum cumulative water deficit (wd), precipitation (pr), temperature maximum (tx), temperature mean (tm), evapotranspiration mean (em), and evapotranspiration mean (et).
# install.packages("geobr")
# install.packages("remotes")
# remotes::install_github("r-tmap/tmap")
library(raster) # Reading, writing, manipulating, analysing and modelling of spatial data
library(sf) # provides simple features access with a geometry list-column
library(tidyverse) # collection of R packages for data science
library(tmap) # library for plotting maps
library(geobr) # library for Brazilian data sets
# setting the working directory
setwd("/data/GY3440/fernando/week5/")
view(geobr::list_geobr())
states <- geobr::read_state()
map1 <- tm_shape(states)+
tm_polygons()
map1
map2 <- tm_shape(states)+
tm_polygons(fill = 'name_state')
map2
legal_amazon <- read_amazon()
biome_amazon <- read_biomes() %>% filter(name_biome == "Amazônia")
map3 <- tm_shape(legal_amazon) +
tm_polygons()+
tm_shape(biome_amazon) +
tm_polygons()
map3
tmap_mode("view")
map4 <- tm_shape(states)+
tm_polygons(fill = 'name_state') +
tm_shape(legal_amazon) +
tm_borders(col = 'green')+
tm_shape(biome_amazon) +
tm_borders(col = 'red')
map4
il <- geobr::read_indigenous_land() %>% st_make_valid() %>% st_intersection(legal_amazon)
pa <- geobr::read_conservation_units() %>% st_make_valid() %>% st_intersection(legal_amazon)
urban <- geobr::read_urban_area() %>% st_intersection(legal_amazon)
map5 <- tm_shape(legal_amazon) +
tm_borders(col = 'green')+
tm_shape(biome_amazon) +
tm_borders(col = 'red') +
tm_shape(il) +
tm_polygons(fill = 'purple', alpha = 0.5) +
tm_shape(pa) +
tm_polygons(fill = 'darkgreen', alpha = 0.5) +
tm_shape(urban) +
tm_polygons(fill = 'grey', alpha = 0.5)
map5
tmap_mode("plot")
map5 +
tm_graticules() +
tm_compass(position = c("left", "top"))
Let’s open the data called prodes_amazonia_raster_2000_2022_v20231109.tif under the folder “data/prodes/prodes_amazonia_raster_2000_2022_v20231109.tif”.
deforestation <- raster::raster("/data/GY3440/fernando/week5/data/prodes/prodes_amazonia_raster_2000_2022_v20231109.tif")
Because this is a TIF (an image) file, we need to add the colour templates to show the deforestation regions. The colours are the below, and the classes are from 2000 to 2020. So, 20 years of deforestation data of the Amazon!
colours <- c("#faef1c",
"#dd3027",
"#e34832",
"#faf82f",
"#e34832",
"#e54e34",
"#e65437",
"#e85a3a",
"#e9603d",
"#eb663f",
"#ec6c42",
"#ee7145",
"#f07748",
"#f17d4a",
"#f3834d",
"#f48950",
"#f68f52",
"#f79555",
"#ffc700",
"#fffebd",
"#feffbf",
"#fafdbe",
"#f7fcbd",
"#f4fbbb",
"#f0f9ba",
"#edf8b9",
"#eaf7b8",
"#e6f5b7",
"#e3f4b6",
"#e0f3b5",
"#dcf1b4",
"#2649e6",
"#2e8113",
"#fc3ad2"
)
brks <- c(-1 ,
0 ,
4 ,
6 ,
7 ,
8 ,
9 ,
10 ,
11 ,
12 ,
13 ,
14 ,
15 ,
16 ,
17 ,
18 ,
19 ,
20 ,
21 ,
22 ,
50 ,
51 ,
52 ,
53 ,
54 ,
55 ,
56 ,
57 ,
58 ,
59 ,
60 ,
61 ,
91 ,
100,
101)
classes <- c("d2000",
"d2004",
"d2006",
"d2007",
"d2008",
"d2009",
"d2010",
"d2011",
"d2012",
"d2013",
"d2014",
"d2015",
"d2016",
"d2017",
"d2018",
"d2019",
"d2020",
"d2021",
"d2022",
"r2010",
"r2011",
"r2012",
"r2013",
"r2014",
"r2015",
"r2016",
"r2017",
"r2018",
"r2019",
"r2020",
"r2021",
"Hydrography",
"Forest",
"Non-Forest")
raster::plot(deforestation,
breaks = brks,
col = colours,
legend = F)
legend("bottom", legend = classes, fill = colours, cex = 0.7, y.intersp = 0.7, ncol = 7, bg = NULL)
Here we are cropping (cutting) the area using geographical coordenates.
e <- raster::extent(-56.5956 , -53.58015 , -5.028585 , -1.727968)
tapajos <- deforestation %>% crop(e)
raster::plot(tapajos, breaks = brks, col = colours, legend = F)
Here we are cropping the area using a shape file data
bb <- st_read("/data/GY3440/fernando/week5/data/Area_corte_Mateus.shp") %>%
st_transform(st_crs(deforestation)) # Transforming vector's projection to match raster's projection
tapajos <- deforestation %>% crop(bb)
raster::plot(tapajos, breaks = brks, col = colours, legend = F)
Both approaches produced the same results.
Loading and plotting the the data from MapBiomas for the year of 2015
Note that the MapBiomas data is more rich. It has more land-use classes such as pasture, and other types of forest. The previou map from Prodes data has only deforestation data - no additional land-use classes.
url <-"/vsicurl/https://storage.googleapis.com/mapbiomas-public/initiatives/brasil/collection_8/lclu/coverage/brasil_coverage_2015.tif"
landuse <- raster(url) %>%
crop(bb)
codes <- readr::read_delim("/data/GY3440/fernando/week5/data/Codigos-da-legenda-colecao-8.csv",
delim = ";", escape_double = FALSE, trim_ws = TRUE)
legend_vars <- codes %>% filter(Class_ID %in% unique(landuse)) %>% arrange(Class_ID)
tm_shape(landuse)+
tm_raster(col.scale = tm_scale_categorical(labels = legend_vars$Description,
values = legend_vars$Color),
col.legend = tm_legend(title = "Land-use for 2015")
)+
tm_layout(legend.outside = TRUE)
Methods: For our burned area estimation, we used 48 Landsat images (Table S2) from Landsat 5, 7, and 8 between the years 2010 and 2016. These images covered an area of 6.48 million ha and included 14 municipalities in central-eastern Amazonia: Aveiro, Barreirinha, Belterra, Itaituba, Juruti, Mojuí dos Campos, Monte Alegre, Nhamundá, Parintins, Placas, Prainha, Rurópolis, Santarém, and Uruará. We performed pixel-by-pixel unsupervised k-means classifications (MacQUEEN, 1967, Drake and Jonathan, 2012) of each Landsat image with six classes and 10 interactions in ERDAS IMAGE v.16 (2016), to classify primary forest (including both undisturbed and disturbed), secondary forest, burn scars (from the 2015- 16 El Niño-mediated fires), deforested areas, bodies of water, and non-forest (figure 1). We used the following as input variables: spectral bands including the visible to the medium infrared, Normalised Difference Vegetation Index (NDVI), Soil Adjusted Vegetation Index (SAVI), Enhanced Vegetation Index (EVI), and the Normalised Burn Ratio 2 (NBR2). Imagery from Landsat 7 and 8 were used in combination with the panchromatic band (Landsat 7 & 8) to improve their spatial resolution. The classified rasters were then imported and vectorised in ArcGIS v.10.2 (ESRI 2014), where a visual inspection of the automatic classification was made to correct any classification errors. Each individual band and all possible combinations in RGB composites were used to identify classifier errors. Following the correction of these errors, we calculated the cumulative area of primary and secondary forest that experienced understorey wildfires during 2015-16 in the Santarém region (figure 1).
bs <- raster("/data/GY3440/fernando/week5/data/Inc_2015_2016_allweeks_500m.tif")
tapajos2 <- projectRaster(from = tapajos, crs = crs(bs), res = 500, method = 'ngb')
raster::plot(tapajos2, breaks = brks, col = colours, legend = F)
plot(bs, add = T, col= "red", legend = F)
tm_shape(landuse)+
tm_raster(col.scale= tm_scale_categorical(labels = legend_vars$Description, values = legend_vars$Color),
col.legend = tm_legend(title = "Land-use for 2015")
)+
tm_shape(bs)+
tm_raster(col.scale = tm_scale_continuous(values = "red"),
col.legend = tm_legend(title = "Week burned"))+
tm_layout(legend.outside = TRUE)
Loading secondary forest
sv <- raster( "/vsicurl/https://storage.googleapis.com/mapbiomas-public/initiatives/brasil/collection_8/secveg-age/sec_veg_age_2016.tif") %>% crop (bb)
tm_shape(sv)+
tm_raster(col.scale= tm_scale_continuous(values = "green",limits = c(1,36)),
col.legend = tm_legend(title = "Age of secundary forest (2016)")
)+
tm_layout(legend.outside = TRUE)
Precipitation data is assemble as the total for each quarter (trimester) for 2015 and 2016, acroding to the table bellow:
| Codename | corresponding months |
|---|---|
| DJF | December, January, February |
| MAM | March, April, May |
| JJA | June, July, August |
| SON | September, October, November |
# loading rivers as auxiliary data
rivers <- st_read("/data/GY3440/fernando/week5/data/Rivers.shp") %>% st_make_valid()
# List of Monthly Rainfall Rasters
monthly = list.files("/data/GY3440/fernando/week5/data/PRECIP_SUM_TRIMESTRAL/", pattern = '.tif$', full.names = T)
# opening all the listed files as stacked raster
precip = stack(monthly) # stacking
tm_shape(precip) +
tm_raster(col.scale = tm_scale_continuous(limits = c(min(getValues(precip))-1,
max(getValues(precip))+1),
n = 7),
col.free = F,
col.legend = tm_legend(title = "Preciptation by trimester (mm)")
) +
tm_shape(rivers)+
tm_fill(fill = "darkblue")+
tm_layout(legend.outside = TRUE)
roads <- st_read("/data/GY3440/fernando/week5/data/Roads_interp_wgs84z21_v2.shp")
raster::plot(tapajos2, breaks = brks, col = colours, legend = F, main = "PRODES Deforestation")
plot(roads$geometry, add = T)
tm_shape(landuse)+
tm_raster(col.scale= tm_scale_categorical(labels = legend_vars$Description, values = legend_vars$Color),
col.legend = tm_legend(title = "Land-use for 2015")
)+
tm_shape(roads)+
tm_lines("black")+
tm_layout(legend.outside = TRUE)
Using cellular space is a method for assessing the landscape. Each
cell works as a sample of the landscape where de variables are
represented. In the loaded dataset as cells, the cells are
the rows and the variables are the columns.
The varibles are described as the following table:
| variable suffix | Description | Category |
|---|---|---|
af_ |
Active Fire (n°/trimester) | Environmental |
ds_ |
Drought Season Length (days/trimester) | Environmental |
wd_ |
Minimum Cumulative Water Deficit (mm/trimester) | Environmental |
pr_ |
Precipitation (mm/trimester) | Environmental |
tx_ |
Temperature Max. (C°) | Environmental |
tm_ |
temperature Mean (C°) | Environmental |
em_ |
evapotranspiration mean (mm) | Environmental |
et_ |
evapotranspiration total (mm) | Environmental |
ag_p |
TerraClass 2014 agriculture proportion | Anthropic |
pa_p |
TerraClass 2014 pasture proportion | Anthropic |
ua_p |
TerraClass 2014 Urban Area proportion | Anthropic |
pf_p |
TerraClass 2014 Primary forest proportion | Anthropic |
sf_p |
TerraClass 2014 secondary forest proportion | Anthropic |
nf_p |
TerraClass 2014 non forest proportion | Anthropic |
wt_p |
TerraClass 2014 water bodies proportion | Anthropic |
pf_edg_l |
TerraClass 2014 primary forest edge length | Anthropic |
sf_edg_l |
TerraClass 2014 secondary forest edge length | Anthropic |
road_lengt |
road length | Anthropic |
CNEFE_n |
CNEFE communities number in the cell | Anthropic |
dt_c |
numbers of DETER forest disturbance warnings | Anthropic |
dt_p |
proportion of DETER forest disturbance warnings | Anthropic |
cr_area_m |
mean of CAR private land areas | Anthropic |
cr_c |
number of CAR private land areas | Anthropic |
cr_p |
proportion of CAR private land areas | Anthropic |
con_p |
proportion of CAR consolidated areas | Anthropic |
af_total |
Active Fire from jun 2015 to aug 2016 | Disturbance |
bs_total |
Burn Scar from jun 2015 to aug 2016 | Disturbance |
The seasonality of the variables are represented in the codename as the following table:
| Codename | corresponding Season |
|---|---|
_djf15_ |
December, January, February 2014-2015 |
_mam15_ |
March, April, May 2015 |
_jja15_ |
June, July, August 2015 |
_son15_ |
September, October, November 2015 |
_djf16_ |
December, January, February 2015-2016 |
_mam16_ |
March, April, May 2016 |
_jja16_ |
June, July, August 2016 |
_son16_ |
September, October, November 2016 |
cells <- st_read("/data/GY3440/fernando/week5/data/cel_10km_v8_ActiveFire.shp")
view(cells)
tm_shape(cells) +
tm_fill(fill = "bs_total",
fill.scale = tm_scale_continuous(values = "red"),
fill.legend = tm_legend(title = "Prop. Burn Scar 2015-2015")) +
tm_shape(rivers)+
tm_fill(fill = "darkblue",
fill.legend = tm_legend(title = "Rivers")) +
tm_shape(roads)+
tm_lines(col.legend = tm_legend(title = "Roads"))+
tm_graticules(lines = F) +
tm_compass(position = tm_pos_out())
tm_shape(cells) +
tm_fill(fill = "af_total",
fill.scale = tm_scale_continuous(values = "red"),
fill.legend = tm_legend(title = "active fire 2015-2015")) +
tm_shape(rivers)+
tm_fill(fill = "darkblue",
fill.legend = tm_legend(title = "Rivers")) +
tm_shape(roads)+
tm_lines(col.legend = tm_legend(title = "Roads"))+
tm_graticules(lines = F) +
tm_compass(position = tm_pos_out())
af <- cells %>% select(contains(c("af_djf15_",
"af_mam15_",
"af_jja15_",
"af_son15_",
"af_djf16_",
"af_mam16_",
"af_jja16_",
"af_son16_"
))) # selecting only the columns of active fire
plot(af, breaks = "fisher",
lwd = 0.5,
key.pos = 4,
pal = hcl.colors(10, "reds"),
max.plot = 23)
wd <- cells %>% select(contains("wd_")) # selecting only the columns of water deficit.
plot(wd, breaks = "fisher",
lwd = 0.5,
key.pos = 4,
pal = hcl.colors(10, "reds"))
Correlation between anthropic factors, active fire
(af_total), and burn scars (bs_total).
cellland <- cells %>% select(contains(c("ag_p", "pa_p","ua_p","nf_p","wt_p", "af_total", "bs_total"))) %>% st_drop_geometry()
cellland[is.na(cellland)] = 0
psych::pairs.panels(cellland, ellipses=F, scale=F, stars=T, gap =0.2, lm=T, cex.cor=1.2, rug=F, hist.col="skyblue3")
celleco <- cells %>% select(contains(c("pf_edg_l", "sf_edg_l", "pf_p", "sf_p","af_total","bs_total"))) %>% st_drop_geometry()
celleco[is.na(celleco)] = 0
psych::pairs.panels(celleco[,], ellipses=F, scale=F, stars=T, gap =0.2, lm=F, cex.cor=1.2, rug=F, hist.col="skyblue3",)
Correlation between active fire (af_djf15_c) and
environmental factors for December 2014, January 2015, and February 2015
season (bs_total).
celldjf15 <- cells %>% select(contains('djf15')) %>% st_drop_geometry()
celldjf15[is.na(celldjf15)] = 0
psych::pairs.panels(celldjf15, ellipses=F, scale=F, stars=T, gap =0.2, lm=T, cex.cor=1.2, rug=F, hist.col="skyblue3")